/*
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
#if !defined(_Cauchy_Wave_FD_)
#define _Cauchy_Wave_FD_
#define NUMBER_OF_PARTICLES 10
#include <src/math/Operators/all.h>
#include <src/algorithm/Evolution_Algorithm/Crank_Nicolson.h>
#include <src/algorithm/Evolution_Algorithm/Escheme.h>
#include <src/math/CauchyProblems/CauchyProblem_traits.h>
#include <src/math/functor/GaussianFunctor.h>


template<class data>
struct CrankNicolson_Wave : public Crank_Nicolson<Wave_t<data>>
{
	OS_STATIC_CHECK(data::Dim == 2); // not sure of the multidimensional behavior. To test before.
	virtual OS_double CFL() {
		return 1. / get_state_variable()->size();
	};
	virtual bool criteria() {
		OS_double breakpoint1 = jmmath::norm22(Op_iterate - Op_S_k1);
		OS_double breakpoint2 = jmmath::norm22(Op_iterate);
		return breakpoint1/(1.+breakpoint2) < 1e-4;
	};
};


//! Second order scheme for the wave equation \f$\partial_t u = \partial_x v\f$, \f$\partial_t v = \partial_x u\f$.
/*!
    This class initialize the Cauchy problem for the one dimensional wave equation
    \f[
        \partial_t w = \partial_x z; \partial_t z = \partial_x w
    \f];
    with \f$x \in [0,1]\f$. This degenerate hyperbolic equation verifies the Energy conservation \f$\partial_t E(u) = 0\f$, with
    \f$E(u) = \int_{[0,1]} w^2+z^2\f$. It is also equivalent to the scalar form \f$\partial_t^2 w - \partial_x^2w=0\f$.

    We illustrate the behavior of the finite difference scheme
    \f[
        \delta_t w^n = \delta_{x} z^{n+1/2}; \delta_t z^n = \delta_{\bar{x}} w^{n+1/2} = 0
    \f]
    This scheme is a second order one. It is entropy (energy) conservative
    for the entropy functional \f$U^n := U(u^n)\f$, with \f$u^n :=
    (w^n,z^n)\f$ and
    \f[
        U(u^n):= |w^n|^2 + |z^n|^2
    \f]
    The boundary conditions are given by <data> : either Dirichlet or Periodic ones.
    The initial data corresponds to a normal distribution centered in 0.5.
*/

template <class data, class Interface = Graphical_Interface>
		class Cauchy_Wave : public CauchyProblem_traits< CrankNicolson_Wave < data >, Interface >
		{
		public :
			virtual ~Cauchy_Wave(){}
			Cauchy_Wave(){
				initial_conditions(NUMBER_OF_PARTICLES);
			};
			virtual smart_ptr<interface_type> Create_Instance() {
				return smart_ptr<interface_type>(new Cauchy_Wave);
			};

			virtual void set_dim(OS_size dim) {
				OS_DYNAMIC_CHECK(dim == 1,"try to initialise a one dimensional Cauchy problem Cauchy_D_x with an higher dimension");
			};
			OS_size get_dim() {return 1;};
			virtual void initial_conditions(OS_size nb_part)
			{
				data result(2,nb_part);
				result[0] = Gaussian_Functor<data::value_type>()(nb_part);
				result[1].assign(0.);
				set_initial_conditions(result);
			};
		};

template <class data>
struct EScheme_Wave
	: public EScheme< Wave_t<data> >
{
	typedef typename EScheme_Wave type;
	typedef  EScheme< Wave_t<data> > Base;
	type() {
		u_.resize(2,0,0);
	}
	virtual void vstar(indice_type n)
	{
		vstar_[n][0] = u_[n][0];
		vstar_[n][1] = u_[n+1][1];
	};
	virtual void gstar(indice_type n) 
	{
//		gstar_[n] = Op_(0,vstar_[n]);
		gstar_[n] = wave_operator(vstar_[n]);
	};
	virtual OS_double Estar(indice_type n) const {
		return jmmath::norm22(u_[n][0])+u_[n][0]*u_[n+1][0];
	};
	virtual OS_double CFL(indice_type n) {
		return .5/(get_state_variable()->size(1));
	};
	virtual bool criteria(indice_type n) {
		OS_double breakpoint1 = jmmath::norm22(u_[n+1] - S_k1);
		OS_double breakpoint2 = jmmath::norm22(u_[n+1]);
		return breakpoint1/(1.+breakpoint2) < 1e-4;
	};
	mutable Wave<data> wave_operator;
};

//! Third order scheme for the wave equation \f$\partial_t u = \partial_x v\f$, \f$\partial_t v = \partial_x u\f$.
/*!
    This class initialize the Cauchy problem for the one dimensional wave equation
    \f[
        \partial_t w = \partial_x z; \partial_t z = \partial_x w
    \f];
    with \f$x \in [0,1]\f$. This degenerate hyperbolic equation verifies the Energy conservation \f$\partial_t E(u) = 0\f$, with
    \f$E(u) = \int w^2+z^2\f$. It is also equivalent to the scalar form \f$\partial_t^2 w - \partial_x^2w=0\f$.

    We illustrate the behavior of the finite difference scheme
    \f[ \label{wave_e}
        \delta_t w^n = \delta_{x} z^{n} \delta_t z^n = \delta_{\bar{x}} w^{n+1} = 0
    \f]
    This scheme is a third order one : we compute that it is equivalent to
    \f[
        \delta_{t\bar{t}} w^n = \delta_{x\bar{x}} w^{n}
    \f]
    It is thus a fully explicit scheme. Its long time time behavior is remarkably stable, 
	when compared to the second order scheme Cauchy_Wave. Ideed, for a CFL constant equal to one, this scheme computes the exact solution.

	The scheme does not comply the Crank Nicolson framework. We use instead the Entropic scheme one :
    First we identify the entropy conservation. We define \f$U^n:=\frac{1}{2}\left(|w^n|^2+z^n\cdot z^{n-1}\right)\f$, and we check that
    our scheme is conservative for this functional :
    \f[
        \delta_t U^n = w^{n+1/2} \cdot \delta_t w^n +z^{n} \cdot \delta_{\hat{t}}z^{n}=w^{n+1/2} \cdot \delta_x z^n +
        z^{n} \cdot \frac{\delta_t z^{n}+\delta_t z^{n-1}}{2} = w^{n+1/2} \cdot \delta_x
        z^n+z^n \cdot \delta_{\bar{x}}w^{n+1/2} = 0
    \f]
	This entropy functional is not positive defined (the scheme is explicit), and thus is stable under the C.F.L. condition
    of hyperbolic kind \f$sup_i\frac{\tau_n}{x_{i+1}-x_i} < C\f$. This entropy conservation defines the entropy variable
    \f$v^{*,n+1/2}\f$ and the flux \f$g(v)\f$
    \f[
        v^{*,n+1/2}=(w^{n+1},z^n);g(v)=(\delta_x z,
        \delta_{\bar{x}}w).
    \f]

    The boundary conditions are given by <data> : either Dirichlet or Periodic ones.
    The initial data corresponds to a normal distribution centered in 0.5.
*/
template <class data, class Interface = Graphical_Interface>
		class Cauchy_Entropic_Wave : public CauchyProblem_traits< EScheme_Wave < data >, Interface >
		{
		public :
			virtual ~Cauchy_Entropic_Wave(){}
			Cauchy_Entropic_Wave(){
				initial_conditions(NUMBER_OF_PARTICLES);
			};
			virtual smart_ptr<interface_type> Create_Instance() {
				return smart_ptr<interface_type>(new Cauchy_Entropic_Wave);
			};

			virtual void set_dim(OS_size dim) {
				OS_DYNAMIC_CHECK(dim == 1,"Cauchy_Entropic_Wave : only dim =1");
			};
			OS_size get_dim() {return 1;};
			virtual void initial_conditions(OS_size nb_part)
			{
				data result(2,nb_part);
				result[0] = Gaussian_Functor<data::value_type>()(nb_part);
				result[1].assign(0.);
				set_initial_conditions(result);
			};
		};


template <class data>
struct Euler_Scheme_Wave
	: public Euler< Wave_t<data> >
{
//! return \f$\tau^n = 0.5 * \delta_i\f$.
	virtual OS_double CFL() {
		return .5/(get_state_variable()->size(1));
	};
};
//! First order scheme for the wave equation \f$\partial_t u = \partial_x v\f$, \f$\partial_t v = \partial_x u\f$.
/*!
    This class initialize the Cauchy problem for the one dimensional wave equation
    \f[
        \partial_t w = \partial_x z; \partial_t z = \partial_x w
    \f];
    with \f$x \in [0,1]\f$. This degenerate hyperbolic equation verifies the Energy conservation \f$\partial_t E(u) = 0\f$, with
    \f$E(u) = \int w^2+z^2\f$. It is also equivalent to the scalar form \f$\partial_t^2 w - \partial_x^2w=0\f$.

    We illustrate the unstable behavior of the first order finite difference scheme
    \f[ 
        \delta_t w^n = \delta_{x} z^{n} \delta_t z^n = \delta_{\bar{x}} w^{n} = 0
    \f]
	CFL is \f$\frac{\tau^n}{\delta_i} = 0.5\f$.
*/
template <class data, class Interface = Graphical_Interface>
		class Cauchy_Euler_Wave : public CauchyProblem_traits< Euler_Scheme_Wave < data >, Interface >
		{
		public :
			virtual ~Cauchy_Euler_Wave(){}
			Cauchy_Euler_Wave(){
				initial_conditions(NUMBER_OF_PARTICLES);
			};
			virtual smart_ptr<interface_type> Create_Instance() {
				return smart_ptr<interface_type>(new Cauchy_Euler_Wave);
			};

			virtual void set_dim(OS_size dim) {
				OS_DYNAMIC_CHECK(dim == 1,"Cauchy_Euler_Wave : only dim =1");
			};
			OS_size get_dim() {return 1;};
			virtual void initial_conditions(OS_size nb_part)
			{
				data result(2,nb_part);
				result[0] = Gaussian_Functor<data::value_type>()(nb_part);
				result[1].assign(0.);
				set_initial_conditions(result);
			};
		};
#undef NUMBER_OF_PARTICLES

#endif

